Exploring black hole superkicks 
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Recent calculations of the recoil velocity in black-hole binary mergers have found kick velocities of 
~ 2500 km/s for equal-mass binaries with anti-aligned initial spins in the orbital plane. In general 
the dynamics of spinning black holes can be extremely complicated and are difficult to analyze 
and understand. In contrast, the "superkick" configuration is an example with a high degree of 
symmetry that also exhibits exciting physics. We exploit the simplicity of this test case to study 
more closely the role of spin in black-hole recoil and find that: the recoil is with good accuracy 
proportional to the difference between the {I = 2,m = ±2) modes of ^'4, the major contribution 
to the recoil occurs within 30M before and after the merger, and that this is after the time at 
which a standard post-Newtonian treatment breaks down. We also discuss consequences of the 
(l = 2,m — ±2) asymmetry in the gravitational wave signal for the angular dependence of the SNR 
and the mismatch of the gravitational wave signals corresponding to the north and south poles. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 98.80.Jk 



I. INTRODUCTION 

More than forty years after Hahn and Lindquist started 
the numerical investigation of colliding black holes 
a series of breakthroughs starting in 2005 0, S [3 has 
turned the quest for stable black-hole inspiral simulations 
into a gold rush. 

A particular focus of the last few months has been the 
so-called recoil or rocket effect due to "beamed" emis- 
sion of gravitational radiation 0, d, 0| • By momentum 
conservation, radiation of energy in a preferred direction 
corresponds to a loss of linear momentum and the black 
hole that results from the merger thus recoils from the 
center-of-mass frame w^ith speeds of up to a few thousand 
km/s. The velocity of this "kick" depends on the con- 
figuration of the system (e.g., the mass ratio and spins) 
and details of the merger dynamics, but not on the total 
mass (velocity is dimensionless in geometric units) . From 
an astrophysical point of view, the recoil effect is par- 
ticularly interesting for massive black holes with masses 
> 10^ M0, which exist at the center of most galaxies 
and may have a substantial impact on the structure and 
formation of their host galaxies. Observational conse- 
qtiences of large recoil have recently been considered in 
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The largest recoil effects have so far been found T^, 
for a particularly simple configuration suggested in 
based on [isj : equal-mass binaries with (initially) anti- 
aligned spins in the orbital plane. Based on numeri- 
cal simulations for different configurations and a post- 
Newtonian approximation [isj . an estimate of 1300 km/s 
had been obtained for this configuration with maximal 
spin [3]. The kicks found in full numerical simulations 
are however even larger, e.g. 2500 km/s ^ or 1800 km/s 
[Til . [T^ for non-maximally spinning black holes. This 
is of the order of 1% of the speed of light, and can be 
larger than the escape velocity of about 2000 km/s from 
giant elliptical galaxies. Extrapolating current numeri- 
cal results for non-maximal spins to maximally spinning 
black holes predicts recoil velocities of up to w 4000 km/s 



[ll| . Smaller but still significant kick velocities have been 
found for several different types of black hole configura- 
tions [H, [H, [13, [H, [H, [la [m, m . Estimations of the 
probabilities to obtain different kick velocities for differ- 
ent mass ratios and high spins were studied in [2^. 

The parameter space of the inspiral of spinning black 
holes is very large, and although its full exploration will 
require numerical methods, analytical understanding and 
approximations will be crucial to render the task eco- 
nomical. The purpose of the present paper is to obtain 
a better understanding of the physics that leads to the 
large kick results recently observed, and in particular to 
compare with post-Newtonian approximations, and see 
where such approximations are accurate, and where they 
(currently) break down. 

We will refer to a configuration similar to that de- 
scribed in [13, [HI, i-e-, two equal-mass black holes with 
spins anti-aligned and in the orbital plane, as a super- 
kick configuration. The superkick configuration exhibits 
"tt symmetry" , i.e. it is invariant under a rotation by an 
angle n about an axis perpendicular to the initial orbital 
plane. It follows from this symmetry that linear momen- 
tum will not be radiated in the x or y directions, but 
only in the z-direction. As a consequence, the center-of- 
mass will remain at (x = 0, ?; = 0), but can move in the 
z-direction. 

The paper is organized as follows. In Sec. [IT] we briefiy 
summarize our numerical methods, and list the simula- 
tions we have performed. Sec. IIIII analyses several as- 
pects of the dynamics of the "superkick" configurations, 
in particular the comparison with post-Newtonian dy- 
namics and various aspects of the (/ = 2, m = ±2) asym- 
metry. Consequences of this asymmetry for the angular 
dependence of the SNR and the mismatch of the gravi- 
tational wave signals, exemplified by the extreme case of 
the north and south poles, are discussed in Sec. IIVI The 
paper concludes with a discussion section and four ap- 
pendices that contain post-Newtonian equations we use 
in this paper, and a number of small results concerning 
the dynamics of moving-puncture simulations. 
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II. NUMERICAL METHODS AND SUMMARY 
OF SIMULATIONS 

In this section we will summarize our numerical meth- 
ods for evolving black-hole binary spacetimes (largely by 
directing the reader to the relevant references) , and spec- 
ify the numerical simulations we performed. The various 
simulations will be motivated more fully later in the pa- 
per; for now we give an overview for later reference. 

We p erformed numerical simulations with the BAM 
[H, [11] and LEAN ^ codes, with modifications dis- 
cussed in Both codes start with black-hole binary 
puncture initial data [13, [13 generated using a pseudo- 
spectral code [1^ , and evolve them with the x- variant of 
the moving-puncture [1, [1] version of the BSSN 
formulation of the 3-1-1 Einstein evolution equations [32l |. 
The gravitational waves emitted by the binary are cal- 
culated from the Ncwman-Penrosc scalar ^1^4, and the 
details of this procedure for BAM and LEAN are given 
in [13| and [1^, respectively. 

The parameters of our simulations are summarized in 
Table [B Each black hole has mass Mi (with mass pa- 
rameter rrii in the puncture data construction [1H|), and 
the total mass is M = Mi + M2. The black holes have 
a coordinate separation of D. In all runs the punctures 
are placed on the y-axis at y = ±D/2 and given mo- 
menta px and spins S = 0.723Mf = 0.2. The spins are 
aligned with the y-direction, except for the runs in the 
"a-series" , which arc characterized by Sy — ±5 cos a and 
Sx = tS sin a. 

The a-scries and P-serics simulations used modifica- 
tions of the MI configuration described in [l^ . This con- 
figuration was chosen because the results showed clean 
fourth-order convergence and high accuracy. We have 
found that the resolution requirements increase signifi- 
cantly for simulations of spinning black holes, and the MI 
configuration, with a small initial separation and there- 
fore short evolution time, provided a convenient starting 
point for our study; these simulations also capture most 
of the important dynamics that we wish to study. 

The a- and P-series simulations were performed with 
the grid setup Xr?=2[6 x 44 : 4 x 88 : 6][88 : 5.82] in the 
notation of [2J], i.e., the six inner boxes had 44^^ points, 
the four outer boxes had 88^ points, the resolution on 
the finest level is M/88, and the resolution at the outer 
boundary is 5.82M. Convergence tests were performed 
for the a = case (which is the same as the MI configu- 
ration in [13]) with inner-box sizes of 40, 44, 48, and cor- 
responding resolutions. Clean fourth-order convergence 
of the linear momentum radiation flux dPz/dt is shown in 
Figure [T] Also shown is convergence in the puncture sep- 
aration, which is not expected to last beyond the merger 
time of about t = 88M since the separation between 
the two punctures inside the common apparent horizon 
quickly approaches zero [13]. 

Further simulations were performed with larger initial 
separation and with quasi-circular orbit parameters 
(calculated according the prescription given in Ap- 
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TABLE L Physical parameters of the simulations performed 
for this paper. 



pendix[^. These arc indicated D6 (for D = 6A/) and 
D8 (for D ^ 8.2M) in the table. The D8 simulation was 
performed using the LEAN code, while all others were 
performed with BAM. The grid setup for the D6 simu- 
lations was the same as for the a- and P-series, and the 
convergence test referred to later used inner box sizes of 
44, 48 and 52 points. The D8 simulation used a grid setup 
Xv^i [2 X 133 : 1 X 155 : 2 X 133 : 3 X 67 : 9] [44 : ff ] , 
where the innermost three levels with 67 points are 
centered around either hole and follow the motion of the 
puncture. 
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FIG. 1: Convergence plots for the puncture separation r and 
the linear momentum radiation in the ^-direction, dP^/dt ob- 
tained for model a = of the a-series. The plots are scaled 
consistent with fourth-order convergence. After merger at 
about t = 88M convergence in the puncture separation is lost 
(as expected). 

Experimentally we have observed that the resolutions 



3 



used in the a- and P-series simulations are not sufficient 
to obtain clean convergence for evolutions of spinning 
black holes orbiting for longer periods of time. It thus 
appears that the good convergence results for these par- 
ticular series are largely due to the close initial separation 
of the black holes, which results in a rather quick merger 
time of about 88M. When the black holes arc placed 
further apart (or even making the seemingly innocuous 
change of choosing quasicircular orbit initial parameters 
for the same separation as the a-series simulations) con- 
vergence is lost before the black holes merge. We expect 
that fourth-order convergence would be obtained if suf- 
ficiently high resolutions were used, but the extra com- 
putational expense was not necessary for the analysis in 
this paper. 

In the D6 simulations we find that the puncture sepa- 
ration and linear momentum radiation flux dPz/dt con- 
verge well for up to 15Af before merger, as shown in 
Figures [2] and [3l Note that since the waves are extracted 
at Rex = 50M, we need to take into account a time lag of 
roughly 50M when comparing times related to puncture 
motion and wave extraction. These simulations will be 
used only for discussions of the qualitative behavior, and 
for analysis at early times, when we are confldent that 
the results are reliable. Similarly the long D8 LEAN sim- 
ulation will only be used for qualitative comparison with 
post-Newtonian results. 



III. ANALYSIS OF SUPERKICK DYNAMICS 

We begin by describing the dynamics of two black 
holes in a superkick configuration with varying degrees 
of simplicity, in order to build up a clearer picture of the 
physics, and to motivate the simulations and analysis we 
have performed. In the simplest picture we draw (Sec- 
tion |IITB|) the spin decouples from the orbital dynamics; 
a more complex picture includes spin precession effects 
fSection llll Cp . and considering the dynamics in full gen- 
eral relativity (GR) in Section IIII Dl allows us to study 
the merger regime, from which most of the kick effect 
originates. The full GR results can then be directly com- 
pared with PN predictions, which we do in Section HITeI 
We discuss the spin of the final black hole in Section llll Fl 



A. Kick velocity and I = 2, m — ±2 symmetry 
breaking 

As noted before, the superkick configuration exhibits 
"tt symmetry" (0 + tt) , thus linear momentum will 
not be radiated in the x or y directions, but radiation of 
linear momentum in the z direction is allowed, and the 
center-of-mass will only move in the z-direction. 

As in nonspinning equal-mass binary simulations, al- 
most all of the energy is radiated in the I = 2,m ~ ±2 
modes: the maximal relative deviation of the energy in 
those modes from the total energies is roughly 2 %, ne- 
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FIG. 2: Convergence of the puncture separation and dP^/dt 
as functions of time for evolutions of model D6. Results are 
scaled for fourth-order convergence. We see that fourth-order 
convergence is lost in the puncture separation at about t = 
115M, which corresponds to roughly t = 165M in quantities 
from waves extracted at R^x ~ 50M, which is about when we 
see a loss of convergence in dPz/dt. Note that we cut the plot 
at t ~ 175 Af when convergence is lost. 



glecting the contribution from the junk radiation. This 
fact, and the symmetry discussed in the proceeding para- 
graph, leads us to expect that we should be able to di- 
rectly relate the kick in the z-dircction to the imbalance 
between the m = 2 and m = —2 modes, i.e., the differ- 
ence in energy that is radiated toward the "north" and 
"south" hemispheres. Using the special relativistic re- 
lation \p\ = E between the momentum and energy of a 
wave packet (traveling at the speed of light) , we expect 
a relation 

Pz=fx {E22 - E2-2) (1) 

for the radiated momentum in the z-direction, where E22 
and £'2-2 ai'c the energies radiated in the Z = 2, m = ±2 
modes, and / is a geometric factor. Here < / < 1 
expresses the fact that the radiation is smeared out in 
solid angle rather than sharply peaked in the direction of 
the poles. Neglecting all modes but I = 2,m = ±2, we 
assume a wave signal in the form ^'4 = kF {t)Y.^'^ {9 , 4>) -f 
\F{t)Y^^^{e,(l)), where k, A are real numbers (k = A in 
the nonspinning equal mass case), F{t) is a complex time 
dependent function, and the ^^2±2 ^'"'^ spin-weighted 
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FIG. 3: Puncture separation and dP^/dt as functions of time 
for the evolutions of model D6. Results from low, medium 
and high resolution simulations are shown. Only the highest 
resolution is shown for dPz/dt. 
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Inserting this ansatz into the expressions for radiated 
en ergy and linear momentum (see e.g. Eqs. (48) and (49) 
in |24l |) we obtain 
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Consequently the value of the geometric factor / can be 
determined as / = 2/3. We find this relationship satisfied 
to very good accuracy in our numerical evolutions, as 
shown in Figure ([3]). 

The relative asymmetry in the energies emitted in the 
I — 2, m — ±2 modes, 2£'22/(-E'22 + -E2-2) (this quantity 
is unity when there is no symmetry breaking) is plotted 
in Fig. O showing a maximal excess of roughly 40%. An 
analytic fit for extraction radius Rext = 50M is 
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FIG. 4: Comparison of the kick velocity in km/s according to 
Eq. (lU, for a range of angles q. Data points for the measured 
kick and the estimate Eq. ([T]) are shown, the points corre- 
sponding to the energy differences are connected. An analyt- 
ical fit, = 2725 cos(176-|-a), to the measured kick is shown 
as a dashed line. Note that Eq. ([1} slightly underestimates 
the kick, which is consistent since it neglects contributions 
from higher order multipoles Z > 2. 



In the extreme case this fit corresponds to E22/E2-2 ~ 
2.4. For this fit the statistical error from 95 % confidence 
in the phase is roughly 20%, and roughly 2% for the 
amplitude of the oscillation. Fits corresponding to the 
extraction radii Rext = 30M and 75M give consistent 
results with Eq. ^ within the statistical error bars. 
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FIG. 5: Excess energy in the I — 2,m — 2 mode, 
2E22/{E22 + E2-2) plotted for extraction radii R^xt = 30 M, 
and 50 M. The curves are the analytical fits for both extrac- 
tion radii, see Eq. ((SJ Clearly, there is no significant depen- 
dence of this ratio on extraction radius. 



B. Simplest assumption: spin decouples from 
black-hole dynamics 

As a first approximation of the dynamics, we may 
imagine that the black holes behave like (force-free) gy- 
roscopes in flat spacetime, and as they orbit each other 
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their individual spin vectors do not change direction. For 
example, if the spins were originally Si = (0, 5, 0) and 
S2 = (0, —S, 0), these vectors would be constant through- 
out the evolution. If we also assume that the spins do not 
noticeably influence the motion in the x-y plane, then 
simulations that start with the same initial separation 
and momenta will display the same dynamics no matter 
how the spins are directed in the orbital plane. This is 
the situation at the first-PN approximation, since spin- 
orbit and spin-spin couplings enter at higher order. 

This picture is surprisingly close to the observed dy- 
namics in numerical simulations. Figure [6] shows the or- 
bital motion in the first six simulations of the a-series, 
each differing only in the initial directions of the spins 
(i.e., Si^2 = S'(t sina, ± cos a, 0) and a = 0...tt in steps 
of 7r/6). The motion shows differences as a is varied 
(shown in the lower panels of the figure), but these differ- 
ences are very small. In contrast, the resulting kick from 
these simulations, shown in Figure IH displays a clear si- 
nusoidal dependence on the angle — the kick varies from 
« -2500 km/s to -1-2500 km/s. (Note that, since these 
values were calculated at a small radiation extraction ra- 
dius Rex = 50M, the values in the plot are systematically 
higher than the correct values by about 10%). Similar 
figures were also shown in [ll[ . 
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FIG. 6: The x- and y-motion of one of the punctures in six 
simulations from the a-series. The dynamics in the x-y plane 
are almost identical for all of the simulations (only one curve is 
actually visible in the upper panels). The small variations in 
the motion (measured with respect to the a = simulation), 
are shown in the lower panels. 

We might conclude from these results that the final 
kick depends only on the initial spin magnitude and di- 
rection of each black hole. One may write an expression 
similar to Eq. (1) in [ll[, which for the superkick case 
reduces to 



Vk = kcos{a — ao), 



(6) 



and determine the constants k and cio, which are ap- 
proximately given by A: ~ 2500 km/s and ao ~ for our 
data. If this simple picture was correct, and the spins re- 
ally did behave as gyroscopes in flat spacetime, then ^ 
would allow us to determine ao = as the spin direction 
that produces the maximum kick. 



C. The effect of spin precession during inspiral 

Needless to say, the real situation is more complicated. 
The picture of the black holes' spins as gyroscopes in flat 
spacetime is valid only at IPN order; at higher post- 
Newtonian orders and in full GR the spin directions 
evolve during the inspiral. We expect that the magni- 
tude of the kick depends on the magnitude and direction 
of the spins when the black holes are close to merger. 
The spin configuration at merger time is a function of 
the initial configuration plus precession effects during the 
evolution. 

The precession effects can be seen in Figure [71 which 
shows {Sx, Sy, Sz) as a function of time for the a = 
simulation described earlier. Here the spins were calcu- 
lated from the black holes' apparent horizons, using the 
coordinate-based method outlined in |33j and also used 
in [2^ . We see that the x- and y-components of the spin 
show noticeable precession during the last orbit, and af- 
ter 80M of evolution the spin vector has rotated by tt/2. 
The z-component shows small oscillations around zero, 
but these are not well resolved with the current accuracy 
of the code. Note however by comparison with Fig. [5] 
that the period of the oscillation is roughly half an or- 
bital period, which is consistent with the post-Newtonian 
equations in Appendix |B] We will further comment on 
the comparison of the spin precession with PN-results in 
Sec. nil El After the formation of a common apparent 
horizon, at about t w 88M, the Sx and Sy quickly drop 
to zero, and Sz jumps to its final value, corresponding 
to the spin of the final black hole, Sz/M'j w 0.7, where 
Mf is the mass of the final black hole. We also com- 
pute the spin of the final black hole from quasi-normal 
ringdown in Sec. [mFl as Sz/M'j « 0.69. Note that this 
value of the final spin is also the value for non-spinning 
equal-mass inspirals, see e.g. f23 |. 

An example for the dependence of the final kick on 
parameters besides a. Figure [8] shows the final kick for 
a = simulations with differing values of the initial mo- 
menta of the black holes (the P-series in TableU]). Chang- 
ing the initial momenta causes the merger time to change, 
and this means that the spins have more or less time to 
evolve, and are therefore in different directions when the 
black holes merge. The initial spin directions are, how- 
ever, the same for all of these simulations. 

This leads us back to Eq. (1) in [llj, which was orig- 
inally written in terms of the spin angles at merger. 
Above we formulated Eq. ^ in terms of the initial an- 
gle, but it would hold equally well if, instead of choosing 
a = a(i = 0) we were to choose some fixed to a-nd use 
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FIG. 7: Evolution of the spins Sx, Sy and Sz of one of the 
black holes over the course of the a = simulation. The 
X- and I/- components of the spin show noticeable precession 
during the last orbit; the spin has rotated by tt/2 after 80M 
of evolution. By contrast, the z-component (lower panel) dis- 
plays small oscillations around zero; these oscillations are not 
well resolved at the current accuracy of the code. Merger oc- 
curs around (at t « 88Af), at which time Sx and Sy drop to 
zero. The spin Sz/M'j jumps to a final value of 0.723 corre- 
sponding to the merged black hole, reflecting the conversion 
of orbital angular momentum into the spin of the final object. 



a = a{t — to) in Eq. ([6]). Only the phase constant ao 
would change. The results shown in Figure[5]suggcst that 
we will also see oscillatory behavior if we make a plot of 
the final kick versus merger time for a series of runs with 
the same initial spin directions. 

In conclusion, given some superkick initial configura- 
tion we need to know both the spin angle a and the time 
until merger in order to predict the magnitude of the final 
kick from initial data. 
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FIG. 8; The final kick as a function of the initial momenta of 
the black holes for simulations with the initial spin angle a = 
0. The differences between the smallest and largest merger 
time is only about 15M; these differences allow a little more 
or less spin precession, and therefore have a strong effect on 
the final recoil, which varies between 2300 and 2700 km/s. 



from a time period of about 60A/. Before that time the 
contribution is negligible, which is even clearer in Fig- 
ure [31 which shows dPz/dt for the quasi-circular model 
D6 which has a longer inspiral phase. We also see in Fig- 
ure [9] that the function dPz/dt obtained for simulations 
with final kick of 2500, and —2500 km/s docs not only 
differ by a mere rescaling. Instead, the curve obtained 
for a = 7r/2 exhibits an additional oscillation. 
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FIG. 9: Plot of dPz/dt as a function of time for the a — 
0,7r/2,7r simulations. Most of the linear momentum is radi- 
ated over a 60M period of time, centered roughly around the 
merger time. 



D. Duration of the recoil 

Simplified models aside, wc know that the final kick is 
due to an integration of dP^/dt over the entire evolution, 
and dPz/dt will be a complicated function of the instanta- 
neous spin directions. A post-Newtonian version of this 
function is given in [l^ and in Appendix [BJ Figure [9] 
shows dPz/dt as a function of time for three simulations 
in the a-series. We make several observations about these 
plots. The main contribution to the final kick originates 



Wc would like to relate dP^/dt to the motion of the 
punctures, but this is not trivial. The radiation is ex- 
tracted at some radius Rex, and plotting this as a func- 
tion of retarded time r — R^x gives only a crude esti- 
mate of what is happening in the vicinity of the black 
holes at any given time. One could try to improve this 
estimate by using instead the luminosity distance (see 
for example |3j|), but we choose to simply look at the 
puncture motion directly. We can calculate the coor- 
dinate acceleration of the punctures in the z-direction, 
az{t) = d^z{t)/dt^. Figure [TUl shows the acceleration of 
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one of the punctures in the a = 0,7r/2,7r simulations. A 
vertical line indicates the time at which a common ap- 
parent horizon forms, and thus gives us an indication of 
how much of the motion is due to effects before merger, 
and how much after merger. 

Although most of the final kick is generated before 
merger. Figure [TOl suggests that dPz/dt is not negligible 
after the merger. Referring back to the plot of dP^/dt in 
Figure [HI the waves from the merger can be estimated to 
reach the radiation extraction sphere between 138M (the 
time when the common apparent horizon forms, 88A/, 
plus the extraction radius, 50M) and 155M (when the 
final black hole's ringdown can clearly be said to have 
begun, by observing that the wave amplitude has a clean 
exponential decay). Whatever time within this range we 
choose to denote as "merger" , it is clear that a significant 
contribution to the final kick arises after that time. 
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FIG. 10: The coordinate acceleration of one of the punctures 
in the z-direction for the a = 0, 7r/2,7r simulations. The ver- 
tical line represents the time of first formation of a common 
apparent horizon. The curves notably differ from those in 
FigureO in particular they show significant extra oscillations, 
while they have roughly the same amplitude The plot essen- 
tially confirms that the gauge dependence of the puncture 
motion in this direction is not well understood (in contrast 
with the orbital motion discussed in Appendix[C]), which pre- 
vents us from drawing quantitative conclusions about the kick 
velocity from the coordinate acceleration. 



E. Comparison with post-Newtonian predictions 

In Figure (TU] we see that there is a significant contri- 
bution to the kick after the black holes merge. Before 
the merger, the main contribution comes from the 30M 
before merger. During this time the black holes arc at a 
separation oi D < 3.5M, and at this separation it is ques- 
tionable whether an accurate post-Newtonian description 
of the radiated linear momentum is meaningful. We will 
now make a comparison between our numerical results 
and the predictions from the 2.5-PN accurate expressions 
in dl. 

Appendices [Bl [C] and |D] summarize the techniques we 
use to compare post- Newtonian and numerical results. 



Briefly, Appendix [B] lists the expressions for the spin 
evolution and radiated linear momenta as found in [13| . 
These expressions were derived in the harmonic gauge. 
Our initial data are instead in the ADMTT gauge (up to 
2PN accuracy [1^), and although it is not obvious how 
well we remain in the ADMTT gauge during evolution 
(but see [s^ for a result that shows excellent agreement 
for larger separations), we would like to see how much the 
results differ between the two gauges. Appendix |D] thus 
gives the expressions necessary to transform the numeri- 
cal quantities, which are assumed to be in the ADMTT 
gauge, to the harmonic gauge. To do this we need to 
calculate the momenta of the punctures as they evolve. 
Appendix [Cl gives 2PN expressions that relate the punc- 
ture's speeds (again assumed to be in the ADMTT gauge) 
to momenta as given in Eq. (jC2p . We see in Appendix [Dl 
that in fact the ADMTT harmonic transformation 
makes little difference to our results over the time when 
the FN approximation is valid. This result may not be 
surprising, but quantifies any confusion that may arise 
when we compare results in the two gauges, and elim- 
inates any major concern that our results may change 
drastically if we were to perform our simulations in the 
harmonic gauge. 

As an aside, these formulas explain the speed at 
which the punctures move in a black-hole binary moving- 
puncture simulation; the puncture speeds and momenta 
are not related by the Newtonian formula p = mv, but 
instead to good accuracy by its 2PN counterpart. 

Before considering the radiated linear momentum, we 
compare the post-Newtonian predictions for the spin evo- 
lution with our numerical results. Figure [Tl] shows the 
evolution of Sx, Sy and Sz for simulation D8, compared 
with the predictions from Eqs. (|B1[) . We see that there 
is very good qualitative agreement in Sx and Sy, even 
close to merger, which occurs at around t = 260Af. The 
z-component does not agree at all well with the 2.5PN 
prediction, but we again note that the puncture motion 
in the z-dircction is much more gauge-dependent than 
the motion in the x-y plane, and it is the positions and 
speeds of the punctures that we use when evaluating the 
right-hand-sides of Eqs. (jBl[) . Furthermore, it is not clear 
how well the numerical determination of the spin based 
on apparent horizons works in this context. Note also 
that the absolute error is very small. The frequency of 
the oscillations of the numerical simulation is rather close 
to the FN result, which is approximately twice the or- 
bital frequency when precession effects are small (cmp. 
Appendix[B|. After merger, a few M after the end of the 
figure, Sz jumps to its final value of around 0.7, and Sx 
and Sy drop to zero. 

In the case of the radiated linear momentum flux 
dPz/dt and the final kick, we know that the assumptions 
underlying the post-Newtonian expressions break down 
when the black holes are very close. Eqs. (jB3[) diverge 
as 1/r^ as the particles' separation r 0, so it is clear 
that a sensible estimate of the kick cannot be made by 
simply integrating this equation. What has been done 
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FIG. 11: One black hole's spin as a function of time for sim- 
ulation D8. Also shown is the 2.5PN prediction for the spin 
evolution, with the puncture dynamics {xi,Vi} used in the 
spin evolution equations (IBlf) . The agreement is very good 
for Sx and Sy, but poor for Sz- This may once again be 
due to the numerical motion in the z-direction being far more 
gauge-dependent than the motion in the x-y plane. At late 
times Sx and Sy vanish, whereas Sz is found to correspond to 
the angular momentum of the final black hole, compare Fig. [7] 
and Sec. HITfI 



in the past (see for example [37|) is to assume a cut-ofF 
separation, and integrate the post-Newtonian expression 
up to that point. We will now show that this approach 
is unlikely to give correct results in the superkick case. 

Figure fT2l shows the function dP^/dt compared to the 
numerical values for the D8 simulation at a retarded time 
t — 54. 5M, chosen to line up the early oscillations in 
dPz/dt, and close to a naive guess of the retarded time for 
the extraction radius Rex = 50M. The post-Newtonian 
values of dPz/dt were calculated as follows. Eqs. (jB3[) 



requires as input the positions, speeds and spins of the 
two black holes. Rather than integrate the full post- 
Newtonian equations of motion, we simply enter the ap- 
propriate quantities from a numerical evolution. This 
allows us to compare, moment by moment, the post- 
Newtonian and numerical predictions of dPz/dt for two 
particles (or black holes) with the {xi,pi, Si} configura- 
tion. 
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FIG. 12: Comparison of the numerical dP^/dt with that pre- 
dicted by Eqs. (|B3p . A time shift of 54. 5M was applied to the 
value calculated from the numerical wave extraction, to ap- 
proximately take into account the wave travel time between 
the punctures and the wave extraction sphere by lining up 
the peak at t ^ 200M. The numerical relativity and 2.5PN 
results agree qualitatively at early times, but diverge quickly 
near merger, probably due to the 1 /r^ term in Eq. (|B3|l . The 
lower plot is a blow-up of the upper plot up to t = 225M. 



In Figure [12] we once again see good qualitative agree- 
ment at early times. At late times the post-Newtonian 
prediction diverges, due to the term in Eqs. (|B3p . 
What is most striking about this plot is that the disagree- 
ment between numerical and post-Newtonian results be- 
comes serious around 50M before merger, which is just 
before the time when the major contribution to the recoil 
begins in Figures [3l and fTOl This suggests that, at least 
in the special case of superkick configurations, if we inte- 
grate the post-Newtonian dPz/dt up to the point where 
its accuracy breaks down, we will grossly underestimate 
the value of the final kick. 

In order to accurately estimate the value of the kick 
analytically, one would need to make a much more so- 
phisticated choice of cutoff separation, and then perhaps 
match to a close-limit analysis. Such a procedure was 
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applied in |38l . |39[ to nonspinning binaries, and may well 
be applicable in the spinning case. One may also be able 
to get good results from a more careful post-Newtonian 
ana,lysis, as was performed (also in the nonspinning case) 
in [40, |41[. We would expect that the superkick case 
would be an extreme and particularly interesting test of 
such methods. 



F. Spin of the final blacl4 liole from ringdown 

In Sec. IIII CI we found that the spin of the final 
black hole as read off from the black hole horizon is 
Jz/M^ w 0.7. Here we will also determine the dimen- 
sionless Kerr spin parameter a = J^/AP from the quasi- 
normal ringdown gravitational wave signal of the slow- 
est decaying spin weighted spheroidal harmonic mode 
(4^ . [4^, which we measure by projecting it onto the 
I = 2,m = ±2 spin weighted spherical harmonics. These 
projected I = 2,m = ±2 waveforms are split into am- 
plitude and phase according to 1^4 = exp(z(p(t)), 
we then perform analytical fits to the waveform for 
170 < t/M < 230, where we see both a clean expo- 
nential decrease of the wave amplitude and a linear in- 
crease of the gravitational wave phase (corresponding to 
a constant frequency). Performing independent fits with 
a linear function for the wave phase and an exponential 
for the amplitude we obtain values for the complex ring- 
down frequency ljqnm- In order to factor out the overall 
mass scale we then perform a lookup of the dimension- 
less quantity Im(wQ7vj\/)/Re(wQ7vj\/) (i-e. essentially the 
inverse quality factor) in a table of QNM frequencies [i^l • 

We will consider in particular data from the a-series, 
see table For the extraction radii Rex = 30M, 

50M, 75M both the Z = 2,m = ±2 results can very 
well be fit with an analytic expression of the form 
ao -I- ai cos(q; + (ySi) + 02 cos(a + '^2), see Fig. (fTS]) . At 
each extraction radius wc get consistent results for the 
amplitudes ao, 0,2 and the phase shifts (^i, Lp2 for the 
m = —2, 2 modes, but we get the opposite sign for 
ai for the m = —2 and m = 2 modes. For ap we 
get (0.6963,0.6891,0.6891) ± 5 x 10"-* (statistical er- 
ror) for extraction radii Rex = (30M, 50M, 75M). For 
the oscillation amplitudes we get consistent results of 
ai = 0.004 ±0.001, 02 = -0.004 ±0.001, with statistical 
errors corresponding to the 95 % confidence interval and 
rounded to one significant digit. We conclude that the 
asymptotic value of the dimensionless Kerr parameter is 
a/M w 0.69, which is consistent with the value 0.7 we ob- 
tained from the black hole horizon. Since the oscillation 
ai cos(q;+(/3i). which has the periodicity of the kick veloc- 
ity, is not consistent between the to = —2, 2 modes, and 
the oscillation 02 cos(q; + (^2) is of the same size, we con- 
clude that both may be non-physical, e.g. they could be 
due to gauge effects in the radiation extraction algorithm 
at finite radius (we suggest an alternative explanation in 
the next paragraph). It is plausible that such problems 
are more serious in the present case of large kicks, than 



when the black hole system does not move with respect 
to the center of gravity. It would be interesting to ana- 
lyze the present case more carefully e.g. along the lines 
discussed in [ist . 

Since the final spin of the black hole is close to the value 
for non-spinning black hole mergers, it appears that the 
individual, anti-aligned spins of the black holes do not 
contribute to the final angular momentum, but rather 
cancel approximately during merger. This is worth not- 
ing since based on the PN analysis and the numerical 
evolutions there is a small oscillating z-component of the 
spin of the individual black holes. At merger time, the z- 
component of the black hole spins is added to the spin of 
the merged black holes due to orbital motion. In princi- 
ple it could happen that the initially small Sz is enlarged 
greatly (as the PN calculation becomes inaccurate), e.g. 
it could happen that the separate spins precess signif- 
icantly towards the z-axis and add significantly to the 
final angular momentum of the black hole. But this does 
not seem to happen, at best there is a small positive or 
negative contribution to the final spin depending on the 
momentary phase of the Sz oscillation during merger, e.g. 
as we see in Fig. [T3l 
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FIG. 13: The plot shows numerical results for J/M^ = 
"22 a2--2 ^pQjjj^g^^ obtained for extraction radius Rex = 50M 
for the final black hole and an analytical fit (solid curve) - 
in the text we conclude that the final Kerr spin parameter 
does not show significant variations (which might be further 
reduced by increased accuracy of the wave extraction). 



IV. EFFECTS OF RECOIL ON SNR AND 
TEMPLATE MATCH 

The gravitational recoil is essentially due to the 
symmetry-breaking between the dominating modes / = 
2, m = ±2. A natural question is how this symmetry 
breaking is reflected in the overlap integrals of the grav- 
itational waveforms. If the symmetry was not broken, 
the gravitational wave signal emitted towards the "north 
pole" would provide the best template also for the south 
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pole. Similarly, for a sequence of waveforms that corre- 
spond to initial data that only differ in spin orientation 
(which we have parameterized by the angle a), we can ask 
how much signal-to-noise ratio (SNR) is lost when try- 
ing to detect the gravitational wave signal corresponding 
to some value of a with a template corresponding to a 
different value of a. 

Answering this question requires accurate waveforms, 
since any mismatch of waveforms can also be due to lack 
of numerical resolution, errors from the finite extraction 
radius (which can be significant, in particular because the 
recoil velocity creates an asymmetry of the geometry of 
the extraction sphere) , and the contribution of the initial 
junk radiation. 

For the data we are considering in this paper, the ini- 
tial junk radiation cannot be separated from the main 
signal in a sufficiently clean fashion, neither can we ob- 
tain accurate error bars on the wave signals of the whole 
Qf-series to really settle the above questions. Nevertheless 
some preliminary results will illustrate the issue. 

We define the correlation function between two time 
series x{t) and y{t) for a time shift r as: 



x{t)y*{t - T)dt, 



(7) 



where a * denotes complex conjugation. Working with 
Fourier transforms 

/oc 
i{f)e^^^f'df. (8) 
-OO 

the correlation function can be written as 

/oc 
iUWU^'^'^'^df (9) 
'OO 

in terms of Fourier transforms x{f),y{f) of the time 
series. The function Rxyir) is thus simply the inverse 
Fourier transform of i(/)y*(/). The value of r for which 
Rxy is maximal determines the time shift required to get 
the maximum correlation between x{t) and y(t). The 
self-correlation R^x is maximal for t = 0: 

/OO 
\i{f)\^df. (10) 
-oc 

More generally we can define the scalar product 

/OO 
iifWim. (11) 
-OO 

The "match" which determines the efficiency of a tem- 
plate y to identify a signal x is defined as 



M — max ■ 



^{x\x){y\y) 



(12) 



to a real time series. We can now evaluate M for sig- 
nals corresponding to different values of the initial spin- 
angle a in our a-series, or for signals corresponding to 
different angles in the sky for a given value of a, say 
one with a large value of the recoil. From symmetry we 
expect that the mismatch 1 — M when comparing the 
signals corresponding to the maximal difference in the 
kick within the a-scries (« 5000 km/s) equals the mis- 
match for the signals that correspond to the north and 
south poles for the maximal recoil case. Indeed we find a 
value of AI w 0.94 ± 0.01 both for the a = case, which 
is close to the maximal recoil and for the maximal mis- 
match case within the a-series. Deviations of ±0.01 here 
come from comparing either the full complex waveform, 
or just hj^ or hx- The uncertainties due to initial junk ra- 
diation, finite extraction radius and numerical error may 
however be larger than 1 %. More accurate data than 
presented here will be required for conclusive error esti- 
mates. Also, a detailed discussion of the dependence of 
the radiation signal on the angle a and the consequences 
for gravitational wave data analysis is beyond the scope 
of the present paper. 

A related question is how much brighter the source 
appears in the direction opposite to the recoil - in which 
more radiation is emitted. For the case of white noise 
we estimate the relative increase in SNR computing the 
ratios of the norm of the strain h for a given inclination 
angle 9 to the strain measured at the south pole {9 = tt) 
by computing: 



SNR(6l) 
SNR(6l = 7 



ih{9)\h{9) 



{h{9 = Tr)\h{9 = tt) 



(13) 



Note that in gravitational wave data analysis the ori- 
entation of a single detector actually reduces the signal 



In Fig. [13] we plot this ratio for the closc-to-extreme case 
member of the a-series a = for h^^ and /i+ — z/ix. 
The excess of signal toward the north pole compared with 
the south pole is roughly 25 %. 



V. DISCUSSION 

We have discussed "superkick" configurations, i.e., two 
equal-mass black holes with spins anti-aligned and in the 
orbital plane, as a simple but extreme "test case" for 
phenomena associated with the large recoil velocities pro- 
duced by spinning black- hole binary systems. The high 
degree of symmetry results in the gravitational wave sig- 
nal being dominated by the Z = 2,to = ±2 spherical 
harmonics, i.e., the recoil is with good accuracy pro- 
portional to the difference of energies radiated into the 
Z = 2,TO = ±2 modes, see Figure 

The asymmetry here is rather strong and in the ex- 
treme case E22/E2-2 ~ 2.4. For gravitational wave de- 
tection the ratio of the amplitude of the strain h is more 
interesting, in the direction opposite to the recoil we find 
an excess of roughly 25 % larger amplitude in the maxi- 
mum recoil case. 
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FIG. 14: The dependence of the expression p3|l is plotted as 
a function of the inchnation angle 9 for h+, hx and h+ — ihx 
for the near extremal member of the a-series a = 0. For 
comparison we also show the curve for the case without spin, 
when h+ is symmetric around 9 = 7r/2. The excess of signal 
toward the north pole compared with the south pole is roughly 
25 %. 

For the large kicks one observes in the "superkick" con- 
figuration, one should certainly worry about the accu- 
racy of the wave extraction. For the present paper we 
have been interested in a qualitative discussion rather 
than very high accuracy, but we point out that a pro- 
cedure to improve the accuracy of wave extraction via 
the Newman-Pcnrose scalar ^1^4 at finite radius has been 
discussed recently in [45j . An overall improvement in the 
accuracy of spinning black-hole binary simulations should 
also be possible by employing higher-order spatial finite 
differencing [4^ and using initial parameters based on 
FN inspiral calculations [47| . 

The main emphasis of this paper has been the compar- 
ison of the dynamics with post-Newtonian predictions. 
We have found that the 2. 5PN- accurate expressions given 
in accurately describe the spin evolution and linear 
momentum radiation up to about 60M before merger. 
After that time the FN estimate of dP^/dt diverges from 
the numerical values. It is also after that time that we 
find the main contribution to the final kick of the merged 
black hole, and this explains why it is difficult to make ac- 
curate predictions of the kick by integrating the FN equa- 
tions up to a cutoff separation. In order to accurately 
analytically model the recoil for superkick configurations 
(and possibly spinning black-hole binary configurations 
in general) we suggest that a more sophisticated post- 
Newtonian treatment would be necessary, or a matching 
of FN methods during the early inspiral with a close- limit 
analysis of the merger and ringdown. It has recently been 
found [23| that a phenomenological formula for the final 



kick, based on the angular dependence of the terms in 
Eq. (|B3p . matches numerical data reasonably well. Hav- 
ing found that the precise form of (|B3[) fails to predict 
the linear momentum radiation in the regime when the 
majority of the linear momentum is radiated, it will be 
interesting to see how well such a phenomenological for- 
mula works for more general configurations, or if a more 
detailed analytic study will suggest a more generally ap- 
plicable formula. 
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APPENDIX A: INITIAL-DATA PARAMETERS 

In the Bowen- York/puncture data that we use, we 
must choose parameters for the masses, separations, 
spins and linear momenta of the two punctures. Most 
of the simulations studied in this paper are based on the 
MI configuration in p^ . for which the momenta were 
chosen without any attempt to have the punctures move 
in quasi-circular orbit (although the eccentricity in the 
resulting evolutions appears to be small). Often, how- 
ever, one wishes to produce quasicircular orbits, and the 
momenta for our non-MI-based simulations are chosen to 
meet that requirement. In this Appendix we describe our 
procedure for calculating those parameters. Note that in 
[47| we have described a procedure to further improve 
the "circularity" of the initial data by using parameters 
obtained from post-Newtonian inspirals. 

In [23] we showed that a 3FN-accurate formula is suf- 
ficient to calculate initial momenta for nonspinning bina- 
ries. In the spinning case we make use of the results of 
Kidder [l3|. They are in harmonic coordinates, but as we 
will see in Appendix IdJ the difference between harmonic 
coordinates and the ADMTT gauge that we expect our 
evolutions to be in are small. Kidder's Eq. (4.7) gives the 
orbital angular momentum of a binary in circular orbit 
as 
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L = MMr)V^LJl + 2f^)-i^ 



X,(L^-S,)(8^ + 777 



3/2 



I 3 

-(5 - 977) - -77x1X2 (Si • S2) - 3(Ljv • Si)(LAr ■ S2) 



i=l,2 



Af2 



3/2 



(Al) 



and the total angular momentum is J = L + S. The 
variables in Eq. (jAip are as follows. The black holes are 
separated by a distance r and have masses M^, the total 
mass is M = Mi + M2, and the mass ratio quantities are 
H = M1M2/M and rj = ji/M . The black holes have spins 
Si, and Xi ~ ■ The quantity Ljv denotes the unit 

vector in the direction of the angular momentum of a 
Newtonian system of nonspinning particles. It need not 
point in the same direction as the full orbital angular 
momentum L of the system, and we may exploit this 
freedom to uniquely find a momentum P that satisfies 



L = 7t(r X P) 



(A2) 



The specific setup of our data is as follows. The punc- 
tures are placed on the y axis and given momenta in the x 
direction. The orbital angular momentum therefore has 
only one component and that points in the z direction. 
In cases where the last term in Eq. (jAip has a compo- 
nent in the x or y directions, we tilt Ljv such that the 
last term is canceled out and L = Lz. Such a case will 
not arise in the situations considered in this paper; the 
last term in Eq. (|Aip will always sum to zero and we can 
simply write px = T^/r. 



APPENDIX B: POST-NEWTONIAN 
TREATMENT OF SPINNING BINARIES 

Consider two particles with masses Mi and M2, spins 
Si and S2, located at positions xi and X2. We define 
X = xi — X2 and v = dx/dt. Expressions for the evolution 
of the spins are given up to 2.5 post-Newtonian order in 
harmonic coordinates by Kidder [l3j . 



S2 X Si 



S. = 1{,L„XS0(2 + |^| 

+3(n-S2)n X Si}, 
S2 - ^ ((L^ X S2) (2 + 1^ ) - Si X S2 



2M2 



-3(n-Si)n X S2}, 



(Bl) 



where the Newtonian orbital angular momentum is given 
by L^v = fJ^i^i- X v). Note that a particular spin supple- 
mentary condition has been chosen. 



The radiated linear momentum is in turn given by 
Newtonian and spin-orbit contributions. 



N = 



8 Sm 2 



105 M \ 7- 



rn 



55^2 - 457^^ 



M' 
12 — 

r 



SO 



50.^-8^ 
r 



{4r{\r X A) - 2w^(n x A) 



15 r5 

^(n X v)[37^(n • A) + 2(v A)]} 



(B2) 



(B3) 



where A = A/(S2/Af2 - Si/Mi) and Sm = Ah - M2. 
Clearly the Newtonian contribution, which was used by 
Fitchett [i^ to provide an early estimate of the recoil 
from the merger of nonspinning binaries, is zero in the 
equal-mass case. 

To evaluate Eqs. (|Bip and (|B3p one needs the particles' 
positions and velocities as a function of time, i.e., one 
needs to solve the post-Newtonian equations of motion. 
Alternatively, we can use the motions of the punctures 
calculated in our numerical simulations. This allows us 
to compare the precession of the spins and the radiated 
linear momentum with that predicted by post-Newtonian 
theory for the same motion. The results of this compar- 
ison are discussed in Section UlI El 

It is instructive to reduce Eqs. (|Bip - (jB3|) to the special 
case of equal mass and 7r-symmetry. Equal masses Mi = 
M2 imply 5m = and hence Pat = 0. Furthermore, 
/i = M/4, T] = 1/4, and A = 2(S2 - Si). 

For TT-symmetry (which implies equal masses), the po- 
sitions of the punctures x^ = {xi, yi, Zi) satisfy 



{x2,y2,Z2) = {-xi,-yi,+zi) 



(B4) 



for all times. This implies for the relative position and 
velocity variables 



x= (2a;i,2yi,0), v = (27;i„ 2i;ij,, 0). 



(B5) 



In words, for Tr-symmetry the punctures can move in uni- 
son in the z-direction, and in general Zi{t) will describe 
an accelerated motion. However, the FN equations (|Bip - 
(jB3|) are expressed in terms of x and v, which describe 
the motion in coordinates in which the center of mass is 
at rest. Both x and v are orthogonal to the z-axis and lie 
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in the z = plane. In particular, v does not have a com- 
ponent in the z-direction, so for general Zi{t) the center of 
mass frame is an accelerated frame. A related statement 
is that in 7r-symmetry the orbital plane remains orthog- 
onal to the z-axis, there is no orbital plane precession, 
and Tjn ~ /^(x x v) is parallel to the z-axis at all times. 
For TT-symmctry the spins can be written as 



Si = C + 



C 



(B6) 



where C and a are the components of the spin parallel and 
orthogonal to the z-axis, respectively. Hence the sum of 
the spins S = Si + S2 points in the z-direction and the 
weighted spin difference A = 2(82 — Si) is orthogonal to 
the z-axis. 



2C, 



-4cr. 



(B7) 



For the initial data we choose C = 0, but in general C, is 
not constant in time. The time derivative of the spins 
can be written as 



7, 



-IN 



S,+3(n.S,)n 



(B8) 

where one part is precession about the z-axis due to the 
orbit-spin term, L^r x Si, but the axis of precession is 
in general not parallel to the z-axis due to the spin-spin 
terms. For 7r-symmetry we obtain 



6 



(n • <T)(n X a) 



|crpsin(2a)z, (B9) 



where z is the unit vector in the z-direction and a is the 
angle between n and a. The z-component of the spins, 
C = S/2, oscillates with sin(2Q;(t)). Since for negligible 
precession a{t) is equal to the orbital phase plus a phase 
shift, we expect two oscillations per orbit, which roughly 
agrees with observation, see Sees. IIII Cl and HTl El For the 
weighted difference of the spins 

A = — 4(T = -[n X v) X cr 

2r^ 

4 

-^(2a + 3|(t| cos(a) n) x C, (BIO) 

which contains precession due to L^r at order l/r^ and 
the spin-spin term at order l/r^. The term Cxcr describes 
a modulation of the precession about the z-axis since C 
oscillates around zero. 

For the radiated linear momentum we note that for 
TT-symmetry the three vectors n, v, and A in (jB3p are 
orthogonal to the z-axis, and hence P50 is parallel to 
the z-axis as it should be. The angle between n and 
V varies slowly over the entire inspiral from about 7r/2 
for circular orbits to a value less than tt for the plunge. 
The angle between the orbital vectors and A oscillates 
with the orbital and precession time scales. Making the 
approximation that n • v « and fi x v « wz, we obtain 



'SO 



7rv(n ■ cr) + Av^(- ■ a) 

\ V 



15 r 



(Bll) 



Even for quasi-circular orbits where in addition we set 
r = there will be radiation of linear momentum in the 
z-direction, which however averages to zero over time. 

Note that in general there are two contributions, one 
proportional to rv and the other to w^, and they are 
offset in phase depending on the angles between the spin 
and the orbital vectors. As the system approaches the 
plunge phase, the rv term should become as important as 
the effects. Note that we have not discussed the Pss 
spin-spin contribution at next PN order, which could also 
be examined for potentially large contributions near the 
plunge, but in numerical simulations of head-on collisions 
the resulting kicks have been found to be small [i^ . 

The PN expressions (|B3j) and the above discussion ap- 
ply in the regime where the post-Newtonian approxima- 
tions are valid. We see in Section [ill El that these expres- 
sions describe the radiation of linear momentum with 
reasonable accuracy up to about 50M before merger. 



APPENDIX C: PN CALCULATION OF 
PUNCTURE MOTION 

In moving-puncture simulations we can readily track 
the motion of the punctures and record their positions 
x{t) and velocities v{t) = ^P{t). We may then be 
tempted to make a Newtonian analogy and guess that 
the puncture's momentum is P = Miv for a black hole 
with mass Mi. However, when we compare this to the 
momentum specified in the initial data, the two values 
differ significantly. For example, evolve two equal-mass 
punctures with initial separation D = 8M, P = 0.14, 
Ml ~ M2 ~ 0.5. From the numerical puncture motion 
we find Miv « 0.075; this value disagrees with P by 
almost a factor of two. 

During a simulation the "punctures" are at an infinite 
proper distance from their black holes' horizons (sol. |5]|. 
and we may worry that correctly physically interpret- 
ing the punctures' motions requires a thorough investi- 
gation of the gauge and geometry of the punctures as 
they evolve. In fact, the punctures' motions can be un- 
derstood from a simple post-Newtonian analysis. 

Up to 2PN order, the Hamiltonian for two point parti- 
cles in the ADMTT gauge and center-of-mass frame has 
been derived in [H, . 

From the Hamiltonian equations of motion. 



dH 
dP' 



(CI) 



where Xi is the separation vector between the two parti- 
cles. At Newtonian order we recover Xi = Pi/{2fi). Up 
to 2PN order we have for circular orbits 



p^r_]_fpHi-3v) 



2n c2 \^ 2^2 

1 /3P'*(1- 577 + 57^2) 



M(3 + 7/) 
R 

MP^{-5 



2O77 + 377^) 
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The more general expressions (removing the assumption 
of circular orbits), and the 3PN terms, will be omitted 
here for brevity. They can be readily calculated from the 
Hamiltonian in [53| . 

As an example, consider the D = 8M quasi-circular 
orbit parameters from the sequence presented in [5^ . 
for which P = O.lllM and the orbital frequency is 
Mil ^ 0.0376. From the orbital frequency we can calcu- 
late that the punctures will move at a speed x = 0.150, 
which is approximately equal to the observed value in a 
simulation. From the momentum, the Newtonian pre- 
diction of the speed is ijv = 0.223, which is far too 
high. The IPN prediction is iipw = 0.126 (now the 
value is too small), and the 2PN and 3PN predictions 
are X2pn = 0.151 and x^pn = 0.149. The 2PN and 3PN 
predictions are both very close to the observed value. 

In addition to providing a pleasing consistency between 
the dynamics observed in moving-puncture simulations 
and that predicted by post-Newtonian theory, this anal- 
ysis is necessary when converting between ADMTT and 
harmonic gauges in Appendix IdJ where we will need to 
invert equations like (|D1[) to estimate the black holes' 
momenta as a function of time from the puncture mo- 
tion. 



APPENDIX D: ADMTT TO HARMONIC 
TRANSFORMATION 



must be determined by the procedure described in Ap- 
pendix o 

When the particles are far apart and moving slowly, the 
coordinates and will not differ much. In Figure [T51 
which shows results from the D6 simulation, we show the 
separation between the punctures in the numerical co- 
ordinates as a function of time, and in the coordinates 
after the transformation (|Dip . The coordinates differ 
by less than 10 % up until about 15A'/ before merger. 
After that time we do not expect the coordinate trans- 
formation (which is accurate up to only 2PN order) to 
be reliable. However, for most of the evolution we see 
that the differences between the two coordinate choices 
are not dramatic. A comparison of numerical and PN 
calculations of dPz/dt (as described in Section fill Ep in 
Figure [11] also shows that the results are similar at early 
times, before the PN result diverges. Note that when the 
puncture motion in ADMTT coordinates is used in the 
(harmonic) PN formula (|B3[) . the curve is closer to the 
numerical result than when we use the puncture motion 
in harmonic coordinates. However, since this agreement 
occurs just before the time when the PN and numerical 
values seriously diverge, we do not take this agreement 
too seriously. 

The main conclusion of this analysis is that the dif- 
ference in results between using ADMTT and harmonic 
dynamical quantities in PN expressions is less than the 
uncertainty inherent in the PN expressions themselves. 
We therefore continue to use the raw numerical data, in 
ADMTT coordinates, for most of the analysis presented 
in this paper. 



Our initial data are, up to 2PN order, in the ADMTT 
gauge [35[. The post-Newtonian expressions for spin 
evolution and linear momentum radiation listed in Ap- 
pendix |B] are in the harmonic gauge. Although it is not 
clear how closely our evolved data adhere to the ADMTT 
gauge, it is nonetheless useful to assume that they remain 
in the ADMTT gauge and transform the results to the 
harmonic gauge and see how different they are. 

A transformation between ADMTT and harmonic co- 
ordinates is provided up to 2PN order by Damour and 
Schafer If are the ADMTT coordinates of the 

i-th particle and are the corresponding harmonic co- 
ordinates, then the transformation for a binary system 
is 



X,: = 



i(n-v,; 



7M, M, 

H 

AR AR 



— V,- V, 

2 * 4 ^ 



The velocities V{a,b} in Eq. 
speeds, but are instead Vi = 



(Dl) 



(jDl|) are not the coordinate 
Pi/rrii, and the momenta pi 
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FIG. 15: Coordinate separation as a function of time for the 
D6 simulation, comparing the numerical data (presumed to 
be in ADMTT coordinates), and the same data transformed 
to harmonic coordinates. The difference is less than 10% up 
until about 15M before merger. It is also clear that after 
this time (when the "harmonic" curve turns upward) the PN 
approximations in the coordinate transformation break down. 
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